function [elements,processInfo] = CalculateMaterialResponse...
    (elements,model,processInfo)
ne = model.countElements();
for i=1:ne
    e = model.getElement(i-1);
    f = e.getFace();
    ue = f.getUHatLocal();
    xi = e.getGaussPoint(0);
    B = zeros(3,6);
    B = e.computeB(f,B,[xi(1) xi(2)]);
    C = elements(i).constitutiveMatrix;
    A = elements(i).A;
    H = elements(i).H;
    Dnn = elements(i).Dnn;
    Dss = elements(i).Dss;
    Dhat = [Dnn 0;0 Dss];
    
    if elements(i).isCrack == 1
        %check crack closure
        %trial value for gamma
        gamma = elements(i).gamma_old;
        %trial value ebar
        e = elements(i).e;
        es = e(2);
        ebar = [0 es]';
        %trial value tbar
        tbar = A*(ue - H*ebar);
        tn = tbar(1);
        ts = tbar(2);
        if tn > 0
            disp(['element ' num2str(i) ' crack opens']);
            %crack opening/re-opening
            elements(i).isOpening = 1;
            elements(i).isSticking = 0;
            elements(i).isSliding = 0;
            %compute etrial
            etrial = (Dhat + gamma*A*H)\(gamma*A*ue);
            etildetrial = computeEquivalentSeparation(elements(i),etrial);
            %compute damage function
            phi = computeDamageFunction(elements(i),etildetrial) - gamma;
            if phi <= 0
                disp(['element ' num2str(i) ' crack unloading']);
                elements(i).isLoading = 0;
                %update separation vector
                elements(i).e = etrial;
                %update gamma
                elements(i).gamma = gamma;
                %update t
                elements(i).t = Dhat*e/gamma;
            else
                disp(['element ' num2str(i) ' crack loading']);
                elements(i).isLoading = 1;
                %update separation vector and gamma
                [e,gamma] = computeGammaOpening(elements(i),ue,etrial);
                %validate e and gamma
                etilde = computeEquivalentSeparation(elements(i),e);
%                 computeDamageFunction(elements(i),etilde)
                assert(abs(gamma-computeDamageFunction(elements(i),etilde))<1e-6,...
                    ['element ' num2str(i) ': error gamma != F(etilde(e))']);
                assert(norm((Dhat+gamma*A*H)*e-gamma*A*ue)<1e-6,...
                    ['element ' num2str(i) ' violate traction equilibrium']);
                elements(i).e = e;
                elements(i).gamma = gamma;
                %update t
                t = Dhat*e/gamma;
                elements(i).t = t;
                %update f
                dFdetilde = computeDamageFunctionDerivative(elements(i),etilde);
                elements(i).f = (dFdetilde*(gamma*t)/Dnn)/etilde;
%                 e
%                 gamma
%                 elements(i).t
%                 elements(i).f
            end
        elseif abs(ts-Dss/gamma*es) + elements(i).mu*tn > 0
            %compressive at crack, residual shear cohesive traction win 
            %friction (i.e sliding mode)
            disp(['element ' num2str(i) ' crack slides']);
            elements(i).isOpening = 0;
            elements(i).isSticking = 0;
            elements(i).isSliding = 1;
            %firstly, assuming damage does not grow
            s = sign(ts - Dss/gamma*es);
            h2 = H(:,2);
            fs = [elements(i).mu s]';
            es = fs'*A*ue/(s*Dss/gamma + fs'*A*h2);
            %check the damage function
            etrial = [0 es]';
            etildetrial = computeEquivalentSeparation(elements(i),etrial);
            gamma_trial = computeDamageFunction(elements(i),etildetrial);
            if gamma_trial < gamma
                disp(['element ' num2str(i) ' crack damage does not grows']);
                %damage actually does not grow
                elements(i).isLoading = 0;
                %update e
                elements(i).e = etrial;
                %update gamma
                elements(i).gamma = gamma; %just to make it clarity
                %update fs
                elements(i).fs = [elements(i).mu s]';
            else
                disp(['element ' num2str(i) ' crack damage grows']);
                %damage grows
                elements(i).isLoading = 1;
                [es,gamma] = computeEsSliding(elements(i),ue,s,es);
                %update s
                elements(i).e = [0 es]';
                %update gamma
                elements(i).gamma = gamma;
                %update fs
                elements(i).fs = [elements(i).mu s]';
                %update f
                etilde = computeEquivalentSeparation(element,elements(i).e);
                elements(i).f = computeDamageFunctionDerivative(element,etilde)...
                    *sqrt(Dnn/Dss);
            end
        else
            %both conditions satisfied: compressive at crack, but crack does
            %not slide (i.e sticking mode)
            disp(['element ' num2str(i) ' crack sticks']);
            elements(i).isOpening = 0;
            elements(i).isSticking = 1;
            elements(i).isSliding = 0;
            %update e, gamma
            elements(i).e = ebar;
            elements(i).gamma = gamma; %damage does not grow in sticking mode
        end
        
        %update strain
        strain = B*(ue - H*elements(i).e);
    else
        %update strain
        strain = B*ue;
    end
    %update stress
    elements(i).stress = C*strain;
end
end

% function [e,gamma] = computeGammaOpening(element,d,e0)
% e = e0;
% etilde = computeEquivalentSeparation(element,e);
% gamma = computeDamageFunction(element,etilde);
% A = element.A;
% H = element.H;
% Dnn = element.Dnn;
% Dss = element.Dss;
% Dhat = [Dnn 0;0 Dss];
% 
% % rhs = A*d - (A*H+Dhat/gamma)*e;
% rhs = gamma*A*d - (Dhat + gamma*A*H)*e;
% 
% conv_criteria = 0;
% while conv_criteria == 0
%     dFdetilde = computeDamageFunctionDerivative(element,etilde);
% 
% %     lhs = A*H+Dhat/gamma - ...
% %         (Dhat*e)*(Dhat*e)'*dFdetilde/Dnn/etilde/gamma^2;
%     
%     f = dFdetilde*(Dhat*e)/Dnn/etilde;
%     lhs = Dhat + gamma*A*H + (A*H*e-A*d)*f';
% 
%     de = lhs\rhs;
%     e = e+de;
%     etilde = computeEquivalentSeparation(element,e);
%     gamma = computeDamageFunction(element,etilde);
%     
% %     rhs = A*d - (A*H+Dhat/gamma)*e;
%     rhs = gamma*A*d - (Dhat + gamma*A*H)*e;
%     
%     if norm(rhs) < 1e-14
%         conv_criteria = 1;
%     end
% end
% end

function [es,gamma] = computeEsSliding(element,d,s,es0)
es = es0;
etilde = computeEquivalentSeparation(element,[0 es]');
gamma = computeDamageFunction(element,etilde);
A = element.A;
H = element.H;
h2 = H(:,2);
Dnn = element.Dnn;
Dss = element.Dss;
fs = [element.mu s]';
rhs = -(fs'*A*d - (s*Dss/gamma + fs'*A*h2)*es);
conv_criteria = 0;
while conv_criteria == 0
    dFdetilde = computeDamageFunctionDerivative(element,[0 es]');
    lhs = -(s*Dss/gamma + fs'*A*h2) + s*Dss*es/gamma^2*...
        dFdetilde*sqrt(Dss/Dnn);
    des = lhs\rhs;
    es = es+des;
    etilde = computeEquivalentSeparation(element,[0 es]');
    gamma = computeDamageFunction(element,etilde);
    rhs = -(fs'*A*d - (s*Dss/gamma + fs'*A*h2)*es);
    if norm(des) < 1e-6
        conv_criteria = 1;
    end
end
end
